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ABSTRACT 

Using numerical techniques we studied the global stability of cooling flows in giant 
elliptical galaxi es. As an init ial equilibrium state we choose the hydrostatic gas re- 
cycling model ( Kritsuk 1996| ) . Non-equilibrium radiative cooling, stellar mass loss, 
heating by type la supcrnovae, distributed mass deposition, and thermal conductivity 
are included. Although the recycling model reproduces the basic X-ray observables, 
it appears to be unstable with respect to the development of inflow or outflow. In 
spherically symmetry the inflows are subject to a central cooling catastrophe, while 
the outflows saturate in a form of a subsonic galactic wind. Two-dimensional axisym- 
metric random velocity perturbations of the equilibrium model trigger the onset of 
a cooling catastrophe, which develops in an essentially non-spherical way. The sim- 
ulations show a patchy pattern of mass deposition and the formation of hollow gas 
jets, which penetrate through the outflow down to the galaxy core. The X-ray ob- 
servables of such a hybrid gas flow mimic those of the equilibrium recycling model, 
but the gas temperature exhibits a central depression. The mass deposition rate M 
consists of two contributions of similar size: (i) a hydrostatic one resembling that of 
the equilibrium model, and (ii) a dynamical one which is related to the jets and is 
more concentrated to the centre. For a model galaxy, like NGC 4472, our 2D simula- 
tions predict M ss 2 M Q yr _1 within the cooling radius for the advanced non- linear 
stage of the instability. We discuss the implications of these results to Ha nebulae and 
star formation in cooling flow galaxies and emphasize the need for high-resolution 3D 
simulations. 

Key words: hydrodynamics - instabilities - cooling flows - galaxies: ISM - X-rays: 
ISM - dark matter 



1 INTRODUCTION 

X-ray observations of normal early-type galaxies with the 
Einstein observatory have revealed an extended diffuse ther- 
mal emission from the hot (0.5-2 keV) compone nt of their 
interstellar medium [ISM] (|Forman ct al. 197§). With lu- 



minosities Lx ~ 10 aa -10 42 erg s" 1 this hot gas contributes 
10 9 — 10 10 Mq to the total galaxy mass in its optical con- 
fines (Forman, Jones & Tucker 1985). The X-ray luminosities 
correlate with the optical (blue) luminosities of early-type 
galaxies, although with a large spre ad of Lx at a fixed Lb 
[Canizares, F abbian o & Trinchieri (1987); Don nelly, Faber 
& O'Connell (]l990|); Eskridge, Fabbiano & Kim (|l995|)]. For 



nearby bright ellipticals the surface brightness distributions 
in X-rays and in the optical are similar as far as the interac- 
tion of the galaxy X-ray emitting gas with the surrounding 
medium is unimportant (Trinchieri, Fabbiano & Canizares 
1986). R OS AT observations revealed correlations of the tem- 
perature of the hot gas with the stellar velocity dispersion 



and with the a bundances for a comp lete sample of 43 ellip- 
tical galaxies ( Davis fc White 1996] ). Galaxies with higher 
velocity dispersions tend to have higher, approximately so- 
lar, abundances and higher diffuse gas temperatures. In all 
cases the gas is substantially (a factor of ~ 2) hotter than 
the kinetic temperature of the luminous stars indicating the 
presence of dark haloes. The analysis of spectral proper- 
ties for a sample of 12 early-type galaxies observed with 
ASCA has revealed systematically lower abundances with a 
mean value of about 0.3 solar and has confirmed the relation- 
s hip between stellar veloc ity dispersion and gas temperature 
(Matsumoto et al. 1997). Gas temperature profiles for el- 
liptical galaxies determined with ROSAT and ASCA have 
been found to be surprisingly uniform for projected radii 
r < 10 r c (r e is the so-called effective or half-light radius). 
They exhibit a positive temperature gradient out to ~ 3 r c 
follow ed by a leveling off or gradu al decrease toward larger 
radii ( Brighenti fc Mathews 1997 ). In addition to hot X-ray 
emitting gas early-type galaxies have been shown to have a 
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cooler comp onent (10 7 — 10 s Mg) detectable at 100 fim (Jura 
is well as warm ionized gas (10 3 — 10 



M,-. 



the central 1-10 kpc) producing optical emission lines (Cald- 


well 19$4; Dcmoulin-Ulrich ct al. 1984 


Phillips et al. 19861 


Buson et al. 1993; Singh et al. 1995 


; Macchetto et al. 



Initially, steady-state cooli ng flow models were invoked 



to inter pret X- r ay observations ( Thomas et al. 1986 
fc Whijel98^ |Sarazin fc White 198$ |Vedder et al. 1988^ 
Sarazin & Ashe 198E). This, in essence, empirical approach 



was based on the simple idea that if the cooling time-scale is 
shorter than the Hubble time near the centre, but still longer 
than the gravitational time-scale, then a very subsonic quasi- 
hydrostatic inflow must take place. However, the global sta- 
bility of these early models has not been proved yet. Later 
time- dependent hydrodynamic modeling has demonstrated 
that, in general, such accretion flows are unstable and suffer 
cooling catastrophes at their centres, which were difficult to 
reconcile with observations ( Vleiksin 1988| ; Murray & Bal- 
bus 19$|). Secular variations of the stellar mass-loss rate 
and supernova heating rate due to stellar evolution in the 
galaxy allowed for a description of the evolution of the hot 
I SM i n ellipticals on the H ubble time-scale [D'Ercole et al. 
(h989|); Ciotti et al. (Il99l|); David, Forman & Jones (1990, 



1991) ]. The resulting models evolved through up to three 
consecutive evolutionary stages: the wind, outflow and in- 
flow phases. Accordingly, most of present-day ellipticals are 
in the outflow phase, while a few of the brightest galaxies 
may have already experienced their transition to the inflow 
regime. Although the problem of the cooling catastrophe was 
not resolved, an important outcome of spherically symmet- 
ric evolutionary modeling was the understanding that the 
use of steady-state cooling flow models to diagnose mass ac- 
cretion rates in dynamical situations can lead to e rroneous 
results dCiotti et al. 199l[ |Murray fc Balbus 199^ ). 

In the absence of direct measurements of the flow veloc- 
ity a systematic study of global hydrodynamic and thermal 
stability properties of the hot ISM in ellipticals can provide 
a selection criterion for realistic flow regimes. Starting from 
the essential physics of cooling flows, we have constructed a 
hydrostatic equilibrium model, which describes a stable in 
situ recycling of gas shed by stars through local thermal in- 
stabilities into a condensate, which can become the material 



for further star formation (Kritsuk 1996). 

In this paper we probe the global stability of the re- 
cycling model in the inflow-outflow context. Section 2 de- 
scribes the basic assumptions of the hot ISM model and of 
the underlying galaxy model. Section 3 gives the details of 
the numerical method used in the simulations. The evolution 
of spherically symmetric perturbations and their stabiliza- 
tion by a conductive heat flux are considered in Section 4. In 
Section 5 we abandon the restrictive assumption of spher- 
ical symmetry and study the development of instabilities 
by mea ns of two - dimensional simulations . This allows for a 



better insight, into the complex physics nf pooling flows The 

discussion and conclusions can be found in Sections 6 and 
7, respectively. 



2.1 Input physics 

It is usually assumed that the source of the hot ISM in 
elliptical galaxies is provided mainly by mass loss in the form 
of winds from red giant stars orbiting in the gravitational 
potential of the stellar system and its dark halo. The gas 
is supposed to be well mixed by 'collisions' of the winds 
originating from stars, which pass by close to each other. 

Supernova (SN) explosions of type la are considered to 
be a major potential heating source for the gas, although 
there is a great uncertainty in the rate of SNe la. SNe can 
also play an important role in the gas enrichment with heavy 
elements, if the ejecta are well mixed with the ISM before 
they cool down considerably. 

We assume that evolutionary processes in the stellar 
system provide a steady source of small-scale non-linear per- 
turbations, which unavoidably become thermally unstable in 
the temperature regime of interest. Heat conduction, which 
is efficient in such a hot tenuous plasma, is able to suppress 
the growth of small- amplitude perturbations in the short 
wavelength limit. However, because of different non- linear 
dependencies of the conduction flux and the radiative cool- 
ing rate on local physical conditions in the gas, finite ampli- 
tude small-scale perturbations will still be able to grow non- 
linear. Hence, a fraction of the gas will quasi-isobarically 
cool down to at least T ~ 10 4 K (becoming > 10 3 times 
denser) and will drop out from the flow, thereby enhancing 
the overall cooling of the gas and removing a fraction of its 
momentum. 

The final fate of the condensing material is still a matter 
of debate. If at some later point star formation occurs in the 
cold condensate, one can expect that feedback heating by 
SNe II and Ib/c will return part of the energy back to the 
hot phase (with some delay). Currently our simplified model 
does not incorporate this feedback heating. 

Thus, we introduce a set of sources of mass, momen- 
tum and energy of the hot ISM, which are distributed over 
the whole galaxy. Their efficiencies depend both on the local 
physical conditions in the gas and on the local parameters 
of the stellar system. It is one of the purposes of our hydro- 
dynamic modeling to follow the evolution of a hot corona 
driven by the competition between sources and sinks. In 
several other evolutionary calculations the rates of processes 
directly related to stellar evolution (e.g., the SN rate) were 
assumed to be variable on a time scale of a few Gyr. How- 
ever, we keep the rates constant and concentrate mainly on 
the short-term evolution of the cooling flows excluding the 
early phases of galaxy formation from our consideration. 

Our hydrodynamic description is valid only on scales 
I > C = max{Z wn , l sn , which are larger than the max- 
imum scale associated with discrete physical sources, such 
as a region blown by a stellar wind (Z wn ), a SN remnant 
(P n ), or a cooling condensation (/*'). The scale C d epends 



on the dis tance from the galactic centre. We refer to ( Math- 



ews 199C|) for a detailed discussion of these issues [see also 
( |Kritsuk 19921 ) for details on thermal instabilities 



Below we briefly define the parameters, which control 
the efficiency of the sources. 



2 THE MODEL 
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2.1.1 Stellar mass loss 



2.1.5 Thermal conduction 



The rate of gas supply to the ISM due to stellar mass loss 
and SNe la is defined as ap*(r), where p* is the stellar mass 
density and a — a*+a sn = const. It incorporates a contribu- 
t ion from 'quiescent' stella r mass loss q t ~ 4.7 x 1 0~ 20 s" 1 
( [Faber fc Gallagher 197^ ; |Sarazin fc White 1987| ) with a 



small addition by supernovae a sn = 3.1 x 10 s 1 (the 
estimate corresponds to r sn - 
and (M/ Lb)* = 8 M Q /L Q ) 



0.55 /i? 5 SNu, M sn = 1.4 M , 
According to van den Bergh 
& Tammann (1991) the SNe la rate, r sn , is uncertain by 



a factor of the order of 1.5, but a significantly lower esti- 
mate for E/SO galaxi es, r B n = 0.15 ± 0.06 h% 6 , was derived 
by Cappellaro et al. ( p.997j ). f\ 



2.1.2 Radiative cooling 

We use a non-equilibrium radiative cooling function 
A(T, Zq) erg cm g s for optically thin hot plasma, de- 
fined in the te mperature range T g [10 4 , 1 s 5 ] K, which we 



compiled from (Sutherland & Dopita 1993) for the case of a 
solar gas metallicity Zq, zero diffuse radiation field, and with 



abundance ratios taken from (Anders & Grevesse 1989) 



2.1.3 Mass deposition 

We define the rate of mass deposition to be a function of 
the local conditions in the hot ISM, p t i = b\{n)p, where 
p is the gas mass density. The spectrum of primeval small- 
scale perturbations in the gas and, in particular, its defor- 
mation by heat conduction are described by a parameter 
b = const € [0, 1]. The Heaviside function x and the linear 
instability growth rate n are defined as 

" n, if n > 0; 
0, otherwise; 



X{n) 



(1) 



d / 5p\ _ 1 /2pA dA\ ap, 



dt 



V T 



dT7 



Here c p is the specific heat at constant pressure. The first 
term in the rhs of equation (^) coincides with Field's insta- 
bility criterion, the second one de scribes stabilization due 
to stellar mass loss (Kritsuk 1992). We assume b = 0.5 for 



the models described in this paper. T he motivation for this 
choice can be found in (Kritsuk 1996). 



2.1.4 Distributed heating 

The rate of heating due to thermalization of stellar winds 
and due to SNe la is assumed to be equal to ap,To, where 
To = (a*T* + a sn T an )/a is the characteristic temperature 
of the heat source. Heating by thermalization of winds from 
stars orbiting in the galactic potential is defined by T* = 

= 6.47 x 10 6 K (we assume <r, = 300 km/s and a mean 
molecular weight p = 0.63; 1Z is the gas constant). The SN 

temperature T sn = ^ = 1.09 x 10 9 K for E sn = 6x 10 50 erg 
and M sn = 1.4 M©. With these values T ~ 7 x 10 7 K, while 
for r sn = 0.15 SNu the source temperature is lower: To ~ 
2.6 x 10 7 K. We assume To = 5 x 10 7 K in our simulations. 



* h 75 = H /75 km s _1 Mpc _1 . 
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The large-scale heat flux q= —rj ■ kVT is defined through a 
reduction factor n < 1, which is treated as a free parameter, 
and the Spitzer ( 1962 ) conductivity coefficient for a fully 
ionized gas 



1.84 x 10 _5 T C 5/2 
K = rnX 6rgSS 



(3) 



where for T > 4.2 x 10 5 K the Coulomb logarithm is In A c 
32.0 + ln[7iT 1/2 (T c /10 7 K)]. 



2.2 'Thermodynamic' equilibrium 

The combination of distributed sources presented above al- 
lows stable steady-state solutions, wh ich were anal yzed in 
the framework of a 'closed-box' model (Kritsuk 1996). These 
solutions describe the ISM in 'thermodynamic' equilibrium, 
i.e. in a state, where the sources and sinks of mass and en- 
ergy are locally balanced. 

An equilibrium temperature of one such state, cal- 
culated for our assumed set of parameter values, is 



T eq (a, b, To, Z) w 1.4 x 10 K. This temperature is close 
to the temperature estimates for hot gaseous coronae based 
upon the X-ray data of nearby giant ellipticals. 

The characteristic time scale, on which the gas is able 
to reach thermodynamic equilibrium from any point in the 
huge attraction region of the phase plane (p,T), 

t s = VtJ^, (4) 

is the geometric mean of the local cooling time t c = 
c2 /[(7 ~ l)/ 9 -^] an d the time scale for stellar mass loss 
t a — (apt/p)" 1 ; c is the isothermal sound speed, and 7 = | 
is the ratio of specific heats. 



2.3 Hydrostatic equilibrium 

For simplicity we assume that the source parameters a, b, To, 
and Z do not depend on radius. This implies that the equi- 
librium temperature T oq is also radius independent. Hence, 
if the gas is in thermodynamic equilibrium, it can hydro- 
statically fill only the potential well of an isothermal sphere. 
It is known, however, that due to dissipation the shapes of 
density profiles of stars and dark matter usually do not co- 
incide in giant galaxies. We therefore choose a potential of a 
two-component isothermal sphere, which incorporates stars 



and dark matter with different velocity dispersions ( Kritsuk 



199f; Kritsuk 1997a). We define the parameter (3 as a ratio 



of these dispersions, 



< 1. 



(5) 



This two-component model is completely defined by four 
parameters. Once, we fix a^, M — ^jT cq (since hot gas and 
dark matter are u sually distributed similarly) and /3=0.5 
[see ( |Kritsuk 1996 ) for details], we have to define only two 
additional quantities, namely, the characteristic radius of the 
gravitating matter distribution 



ro = 



(6) 
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and the ratio of dark matter density to stellar density at the 
centre 



S — Pdm,o/p* 



(7) 



We set S = 10 -1 ' 5 , since the stellar component dominates by 
mass in the central regions of giant ellipticals. Such a choice 
of /3 a nd 8, which con trol the shape of the density distribu- 
tions (Kritsuk 1997b), corresponds to a galaxy, whose sur- 



face brightness profile follows the de Vaucouleurs law, and 
which contains comparable amounts of mass in luminous 
and dark matter within its effective radius. Finally, in order 
to fix the scaling of physical values, we set the characteristic 
radius ro = 0.16 kpc. This particular choice of parameters 
allows us to reproduce with reasonable accuracy the shape 
of the surface brightness profile for NGC 4472 (M49), the 
optically most luminous E2 g alaxy in the Virgo clus ter, as- 
suming a distance of 17 Mpc ( rrinchieri et al. 1986 ). 

Numerical integration of the Poisson and Jeans equa- 
tions for the given set of parameters provides us with a 
giant spheroidal model galaxy with a total mass in stars 
M* ss 8.5 x 10 11 Mq and a stellar velocity dispersion 
(T* « 304 km s _1 . 

The hydrostatic isothermal gas distribution in its grav- 
itational potential follows the so-called /3-law: p cq oc p^ with 
f3 — 0.5. This implies that, as it is observed, the shapes of 
the optical and the X-ray brightness profiles are identical 
[cf. ( [rrinchieri et al. 1986) )] and T cq = 2T* [cf. ( |Davis & 
White ]L996| )]. The self-gravity of the ISM is negligible and 
is not taken into account. 

The central equilibrium gas density is p eq ,o / (pmu) ~ 
0.19 cm -3 and the time scale associated with the sources at 
the centre is t a ,o ~ 4.9 x 10 7 yr. The gravitational time scale 



to 



1 _ ro 



2.1 x 10 b yr 



(8) 



is much shorter than t Sl Q. The cooling radius, fixed at the 
distance where t c — 10 Gyr, is r coo i = 36 kpc. 



2.4 Basic equations 

We probe the global stability of the equilibrium hydrostatic 
gas configuration with respect to a variety of perturbations 
by means of numerical solutions of the time-dependent hy- 
drodynamic equations for the mass, momentum and energy 
balance in the hot ISM: 

dp 
at 



+ V • (pv) = ap* - p ti , 



(9) 
(10) 



+ V (p« 2 + p) = pV0 - ptiv, 

^+V.[v(E + p)] = 

otp*e» - p\i(E + p)/p - p 2 A + pv ■ V0 - V • q. (11) 

Here v is the velocity, p — (7 — l)ep is the pressure, e is the 
specific internal energy, e* = c v Tb, and the energy density 
E = p(e + v 2 /2). 



2.5 Initial conditions 

The undisturbed state is considered to be a spherically sym- 
metric gas configuration in hydrostatic and thermodynamic 



equilibrium: 

P| t = = Peq(r), 



'R 



P\ t = = Pcq(r) = — pcqTeq, 

vl_ n = 0. 



(12) 
(13) 
(14) 



For one-dimensional spherically symmetric simulations we 
introduce a global density perturbation with an amplitude 
e = const: 



p| t=0 = (l + e)p eq (r). 



(15) 



For the 2D axisymmetric simulations, which were performed 
in spherical coordinates (r, 1?, ip), we used either a single- 
wave isothermal density perturbation of the form 

p| 4=0 = [1 + e (cos 2 - 0.5)] p oq (r), (16) 

or random large-scale velocity field perturbations, v = 
(u,v), (where u and v are the radial and the tangential 
velocity component, resp ecti vely) which are defined on the 
discrete grid (see section 3.3) as follows: 



u,v\ t=0 = <Ju,«fV,0), 



8 a (rk,$i) 



1 3 3 

i=i j=i 
2tt 



( 2tt . 2tt 

An cos — !hm — 11+ 



(17) 



(18) 



2tt 



B° cos — ik cos — jl + 



2tt 



2tt 



Cy sin — ik sin — jl + 



2tt 



2tt 



D"j sin — ik cos — jl J , 



with a £ {u, v}. The computational volume is discretized 
into n x n cells, which are labeled by the pair of indices k, I. 
{A, B, C, D}tj are quasi-random numbers from the interval 
[—1, 1]. The amplitude of the velocity perturbations is nor- 
malized by a factor Q, so that i(maxfc,; 5 a — mint,; 5 a ) = e.f] 
We varied the amplitude e in different experiments in 
the range from 0.2% to 30%. 



2.6 Boundary conditions 

We imposed the following conditions on the hydrodynamic 
variables at the inner spherical boundary at ri w 0.12 kpc: 



u{n) =0, — 



«(n) = 0, 1 r- 



0, 



0. 



dT 

Or 

dp 
Or 



= 0, 



(19) 
(20) 
(21) 
(22) 



At the outer spherical boundary at r w 151 kpc, we used 
slightly different conditions: 



t Velocities are given in units of the adiabatic sound speed 
c(T ) ft! 812 km s- 1 . 
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u(r a ) = 0, 
v(r ) = 0, 



du 
dr 

du 

dr 



= 0, 
0. 



dT 

T(r ) = T cq , — 



p(r D ) =p eq (r ), 



9p 



(23) 
(24) 
(25) 
(26) 



In the 2D axisymmetric calculations we solved the flow equa- 
tions in a 90 degree sector centered at the equatorial plane 
and imposed 'periodic' boundary conditions in angular di- 
rection at $i r = 5 ± ? : 



/- q \ da 



a(#i 

where a G {p, u, 



9o 



(27) 



3 NUMERICAL METHOD 

The numerical solutions discussed in this paper were ob- 



tained using the Pi ecewise Parabolic Method (PPM) ( Colella 
& Woodward 1984 ) - a second order extension of Godunov's 
finite difference scheme, which uses a nonlinear Riemann 
solver in the hydrodynamic calculations. It was orig inally 



described in (Godunov 1959; Godunov et al. 1961) and 



first a dapted for high order difference schemes by van Leer 
(1979). We use a direct single-step explicit Eulerian formu- 
lation of the method for the equations written in conserva- 
tion form. We solve the equations in spherical coordinates 
for both spherically symmetric and axisymmetric problems. 
The coordinate singularities at the centre and along the sym- 
metry axis are isolated by a proper choice of the computa- 
tional volume and the corresponding boundary conditions 
(see above). 

There are slight differences in our PPM implementa- 
tion, as com pared to the original prescription of Colella & 
Woodward ( 1984 ) . These are stipulated by the presence of 
source terms, which play an essential role in our modeling. 



3.1 Treatment of the source terms 

In order to keep the numerical scheme consiste nt, the ex- 
pressions for 0° and [see equation (3.7) of (Colella & 
Woodward 1984)] have been re-derived to incorporate the 
terms corresponding to the volume sources of mass, momen- 
tum and energy, due to radiative cooling, stellar mass loss, 
SNe heating and local thermal instabilities. This gives rise 
to a modification of the procedure for the calculation of the 
effective left and right states for the Riemann problem. The 
effect of the source terms on is accounted for to 

sec ond order acc uracy. The detailed formulae can be found 
in ( |Kritsuk 1998 



The presence of sources, represented by non-linear func- 
tions of the hydrodynamic variables, considerably modifies 
the behaviour of the gas flow. In particular, when the ef- 
fects of local thermal instabilities are taken into account, 



t We keep here the notation of Colella & Woodward ( 1984 ) 
© 1997 RAS, MNRAS 000, [|-|l| 



the elastic properties of the gas change dramatically. In sim- 
ple terms, slow gas compression by an external force easily 
results in condensation of a significant fraction of the gas. 
The self-accelerating contraction does not meet any consid- 
erable resistance by pressure forces, as long as it is slow 
enough. When simulating such gas flows we observed the 
formation of 'condensation fronts' - thin sheet-like quasi- 
isobaric structures, characterized by an enhanced density 
and negative velocity divergence. These fronts, being the 
major sites of condensation, are able to move in the direc- 
tion opposite to large-scale pressure gradients and coalesce 
with each other upon collisions. In order to treat this kind 
of solutions properly with our numerical scheme, we had to 
introduce additional dissipation in the vicinity of the fronts. 
The dissipation algorithm allows us to keep the peak density 
values finite, limiting the run-away cooling and condensation 
of the gas in cells located at the front. In fact, we added an 
explicit diffusive flux to the numerical fluxes in the mass 
and momentum conservation equations following the simple 
prescription of Colella & Woodward (1984). Further detail s 
of the dissipation algorithm can be found in (Kritsuk 1998). 
While we apply the simple flattening around s hock fronts 
[see the Appendix in ( Colella fc Woodward 1984 )] in our ID 
calculations, we have switched it off in our 2D simulations. 
The sources make the method implicit, since the values 



ofp" +1 and p 



n+l, 



_1 depend on E n+1 through terms, which 



include temperature dependent radiative cooling A and mass 
deposition pu. To avoid an iterative solution method we ap- 
proximated the unknown values at t — t n+1 by those from 
the previous time level at t — t" in the source terms p" ; +1 
and A n+1 , still using the correct treatment elsewhere. 



3.2 Introducing the heat flux 

In order to estimate the stabilizing effects of thermal con- 
ductivity on the isothermal equilibrium solution we include 
a V • q term in the final conservative difference step. To keep 
the scheme explicit, we estimate the temperature values at 
zone interfaces using the computed left and right state val- 
ues for pressure and density. Linear approximations of the 
temperature gradients at the interfaces are computed using 
the temperature differences in adjacent zones. An adequate 
reduction of the time step was implemented, too. Such a 
simplified treatment for the heat conduction flux still allows 
us to keep the accuracy of the computations sufficiently high 
in a broad vicinity of the equilibrium solution. 



3.3 The computational grid 

The computational grid is organized in such a way that it 
allows us to describe the physically discrete sources as con- 
tinuously distributed in space. In our particular model the 
largest scale associated with an individual source is the scale 
for local thermal instability (l t[ > l sn > P w ). This scale de- 
pends on the degree of suppression of heat conduction by 
magnetic fields and on the spectrum of the primeval per- 
turbations, which are supplied by ongoing evolutionary pro- 
cesses in the galaxy. We introduce a correction factor rj, 
which defines the scale for growing local thermal instabili- 
ties l u ~ 0.3r/Tf /ri-i kpc. With r\ w 0.1, in the core region 
of the gas distribution I ~ 0.03 kpc becomes smaller than 
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o.i 1 10 io 2 o.i l 10 io 2 

r (kpc) r (kpc) 

Figure 1. Inflow solution t = 2.4 X IO 7 yr after the start of the simulation for an e = 10 % global spherically symmetric density 
perturbation (solid lines). Dotted lines indicate the unperturbed hydrostatic equilibrium solution. The plots show density, temperature 
and radial velocity of the gas (left column) together with the surface brightness in X-rays, J x , the mass deposition rate and the efficiency 
of mass dropout due to thermal instabilities (right column) as functions of the radius. The surface brightness is given in arbitrary units, 
see text for details. 



the minimum scale on which supernovae contribute as a con- 
tinuously dist ributed sourc e, l sn ~ 0.06 kpc for the assumed 
SN rate [see ( Kritsuk 1992 ) and references therein for more 



details]. Thus, if the grid cells are larger than w 10 kpc in 
the core region, the sources can be treated as being contin- 
uously distributed. 

The grid is organized as follows. The inner k radial cells 
are equidistant: Tj = rj-i + A r for j — I,. . . ,k— 1, starting 
from n « 0.12 kpc, where A r « 0.020 kpc, and k = 20. 
Hence, the core region of the initial density distribution 
(King's core radius is r c ~ 3ro) is covered by a uniform 
radial grid. At larger radii an exponential grid spacing is 



used: r$ — rj_i A r exp (6.2- 



for j = k, 



1. This 



spacing guarantees a smooth transition from a uniformly 



to a non-uniformly spaced grid at r ~ 0.40 kpc. The total 
number of radial grid points is m = 128. 

In the 2D simulations we used the same r-spacing and 
an equidistant angular grid of 128 zones with = 7r/2/128. 



4 EVOLUTION OF SPHERICALLY 
SYMMETRIC PERTURBATIONS 

4.1 Case of suppressed heat conduction 

First, we discuss the global stability of the hydrostatic equi- 
librium gas configuration assuming that heat conduction is 
completely suppressed on large scales. 

Our simulations clearly demonstrate that in the absence 
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of efficient heat conduction the equilibrium is unstable. For 
example, in response to a 10% positive density perturbation 
of the form @ an inflow develops on a central cooling time- 
scale t c ,o ft* 0.5£ s ,o (see Fig. [I]), whereas the same negative 
density perturbation initiates an outflow (see Fig. ^|). 

With our numerical scheme we are not able to describe 
the hydrostatic equilibrium of the initial gas configuration 
perfectly. Consequently, the 'unperturbed' initial model gen- 
erates some wave motions, which help the gas to settle in 
the potential well. This provides some initial heating, i.e. 
the 'unperturbed' model generates an outflow. The mini- 
mum amplitude which is necessary to overcome the initial 
noise inherent in the numerical approximation of the equilib- 
rium model is e w 0.0045 for spherically symmetric density 
perturbations 



4- 1.1 Inflow solutions 

A positive gas density perturbation first breaks the balance 
between cooling and heating in the core region of the gas dis- 
tribution. Gravity immediately starts to accelerate the gas 
inflow. As a result, a dense nearly isothermal gas core with 
T « 10 6 K forms, which is condensationally unstable. The 
further development of instability depends on the detailed 
structure of the inflow velocity profile. The condensation oc- 
curs first at a radius, where V ■ v has a local minimum. In 
the example shown in Fig. |l| this happens nearly simultane- 
ously at the inner boundary and at r w 0.5 kpc. Two sharp 
minima in the temperature profile (with T ~ 10 4 K) and 
depressions in density distribution due to the high conden- 
sation efficiency can be clearly seen. For r\ € [0, 0.1] the mass 
loss rate 

IT 

M(r) = 4tt / r 2 p ti Ar (28) 
Jo 
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at this transient stage is 2 — 4 Mq yr 1 within r coo i- Most 
of the condensate is deposited inside a radius of 1 — 2kpc. 
The surface brightness ^| in the core region grows by a fac- 
tor of 10 over its equilibrium value. We had to terminate 
the calculation at t — 0.49 t s ,o ~ 2.4 x 10 7 yr because the 
iteration within Riemann solver no longer converged when 
the temperature drops below 10 4 K. 

Such accelerated, non-linear unstable inflow regimes can 
obviously exist only in configurations with perfectly sup- 
ported spherical symmetry. This is an unlikely situation for 
real galaxies, and hence points towards the need for multi- 
dimensional simulations. 

The existence of an exact hydrostatic equilibrium solu- 
tion allows us to define the boundary conditions properly. 
Conditions at the inner boundary are of particular impor- 
tance for the case of inflow solutions. Free boundary condi- 
tions (i.e. zero spatial derivatives of all dependent variables) 
are inconsistent with hydrostatic initial conditions, because 
the gas is forced to flow towards the centre, having no hy- 
drostatic support at the inner boundary. Our choice [see 
equations (^)-(^)] seems to be optimal for this particular 
problem. In case of a fixed temperature at its equilibrium 
value, T | = T cq , spurious oscillations are generated near 
the inner boundary which accelerate the development of the 
central cooling catastrophe considerably. 

A kind of settling 'cooling flow' solution was found by 
David et al. (1990) in time-dependent simulations of hot gas 
flows in elliptical galaxies. They ignored thermal conductiv- 
ity and included both, sinks and sources for the gas mass, 
although with a slightly different parameterization. Appar- 
ently, the existence of these settling solutions is mainly due 
to their specific choice of the inner boundary conditions. 
They use a King profile for the stellar mass distribution and 
keep the derivatives of temperature and density equal to 
zero at the inner boundary, located in the core region, al- 
lowing for a free gas inflow through the boundary towards 
the galactic centre. This simply cuts out the central region, 
where the cooling catastrophe would have developed. 



4-1.2 Outflow solutions 

The slow outflow solution describes the formation of a hot, 
subsonically expanding bubble filling the central r > 10 kpc 
of the galaxy (see Fig. |^). The gas density gradually drops 
below the equilibrium values in the region covered by the 
outflow, while its temperature rises up to T w 4 x 10 7 K. 
Gas deposition due to thermal instabilities in the overheated 
and rarefied plasma becomes inefficient, and the gas shed 



S Here we conventionally call 'surface brightness' the emission 
measure of the unit column along the line of sight. The X- 
ray surface brightness includes comparable contributions from 
both p 2 A and p t ;e, so that the combined spectral volume emis- 
sivity of the cooling hot gas and lo cal unstable c ondensates is 



Table 1. Stabilizing effect of the conductive heat flux" 



2 (r)A„(T) + p u (r)c v T u (T) (Kritsuk 199^). This can be 



rewritten in the form c v (r) = pr(r)fi(y, T) — p\^i(y,T), where 
/l 2 are functions of temperature and frequency. In order to get 
an idea of the shape of the surface brightness profile, which can 

--. One has to be 



be observed, we define h = — f ° 

' A r Jr 



,0^ 



inflow: t cc /*s,0 



outflow: M n 



[km s" 1 ] 



±e n = 

0.005 1.59 

0.01 1.08 

0.1 0.49 



n = 0.3 
oo 
1.76 
0.59 



r; = l 
oo 
oo 
1.33 



n = 
196 
197 
241 



77 = 0.3 n = l 

20.5 2.3 

21.3 2.8 

32.5 10.6 



a The time-scale t cc for the development of the cooling catastro- 
phe is given in units of the characteristic time for sources at the 
centre, t Si o- Infinity signs indicate that the cooling catastrophe 
does not develop. In this case the solutions either show an oscil- 
latory behaviour around the hydrostatic one or evolve into a very 
subsonic outflow. The values of a ma i give the maximum veloc- 
ities reached by the flow during the evolution until t = 10 t s ,o- 
The initial perturbation amplitude e is positive in case of inflow 
solutions and negative for outflows. 



by stars starts to blow out from the centre. As a result, a 
characteristic feature of cooling flows - a strong central peak 
of the X-ray surface brightness - is being smeared out and 
the X-ray luminosity of the galaxy decreases dramatically. 
We terminated the calculation at t ~ 4i Sj o, when the leading 
wave front of the outflow region reaches the outer boundary 
at r Q = 151 kpc. 

Spherically symmetric outflows develop in a very stable, 
gradual way and, obviously, can serve as an explanation for 
the low X-ray-to-optical luminosity ratios Lx/Lb of some 
ellipticals. In our case they originate in a natural way from 
an equilibrium gas configuration and do not require any co- 
ordinated tempora l changes of the S N rate and/or stellar 
mass loss rate, cf. (Ciotti et al. 1991). 



4.2 



Stabilizing the equilibrium by thermal 
conductivity 



cautious, however, recovering the surface brightness profiles from 
the emission measure in essentially non-isothcrmal flows. 



In order to check the effects of the conductive heat flux, we 
made a number of runs, varying the reduction factor r\ and 
the perturbation amplitude e in a wide range of values. 

The data given in Table 1 show that thermal conduc- 
tivity is more effective in damping the non-linear develop- 
ment of negative density perturbations than those of posi- 
tive ones. It efficiently suppresses the development of out- 
flows and keeps the gas velocities quite low even for quite 
large initial perturbations (see Fig. [|). It also affects the 
development of positive density perturbations by delaying 
the onset of the cooling catastrophe and by creating a neg- 
ative temperature gradient in the core region which shifts 
the onset of condensation towards the inner boundary (see 
Fig. We find, however, that a heat flux is able to pre- 
vent catastrophic cooling only, if it is not suppressed and 
if the initial perturbations are sufficiently small. Solutions 
obtained in runs with rj = 1 and e = 0.5% or 1% exhibit a 
stable quasi-hydrostatic behaviour at t « 10 t 3 ,o with devi- 
ations from the equilibrium temperature of the order of 1% 
and inflow/outflow velocities u < 1 km s _1 . 

We further find that even a considerably suppressed 
heat flux damps the small-scale numerical noise, which one 
can infer from the temperature panel in Fig. This noise is 
amplified by a physical instability similar to that occurring 
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o.i l 10 10 s 

r (kpc) 

Figure 3. Example of a spherically symmetric inflow solution at t ■■ 
r\ = \. See also Fig. 1. 




6 X 10 yr. The conductive heat mix is taken into account with 



in Field's instability. If r\ = 0, the noise saturates at a level 
of the order of 1%. 

The results of these simulations confirm the earlier qual- 
itative finding of Meiksin (1988) , who showed that heat 
conduction delays the development of the cooling core, but 
does not prev ent it (see Table 1). However, in contrast to 
(Meiksin 198c), we did not find any solution of the 'cooling 
flow' type with a considerable inflow velocity, which is rea- 
sonably stable. We attribute thi s difference m ainly to stellar 
mass loss, which was ignored in (Meiksin 1988). It provides a 
steady mass supply and simultaneously stabilizes local ther- 
mal instabilities by reducing the sink terms [see equation 
(g)]. Both effects amplify the cooling and destabilize the 
flow near the centre. Since central galaxies do always reside 
in cooling flows, stable inflow solutions of this type seem to 
be unrealistic. 

Instead, our results demonstrate that steady-state slow 



outflow solutions of the type shown in Fig. [| should be com- 
mon in ellipticals, since they are more stable than any of 
the inflow solutions. Their observable properties are close to 
those of the equilibrium model, while the mass deposition 
rate is lower, M(r coo i) ~ 0.4 Mq yr -1 within the cooling 
radius. At the same time, the slow outflow solutions do not 
show a central temperature depression. 



5 AXISYMMETRIC PERTURBATIONS 

If a spherically symmetric hydrostatic equilibrium gas con- 
figuration is unstable it evolves either into inflow or into 
outflow. If one relaxes the assumption of spherical symme- 
try and allows for 2D perturbations and flow pattern, one 
can get hybrid solutions which show both types of unsta- 
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Figure 4. Example of a steady-state slow outflow solution at t = 4.9 X 10 s yr. The conductive heat flux is taken into account with jj = 1. 
See also Fig. 1. 



ble behaviour in one model. Our axisymmetric simulations 
demonstrate that this is indeed the case. 

We simulated the evolution of axisymmetric perturba- 
tions of various kinds in a 'toroida l' vo l um e centered around 
the equatorial plane (see sections 2.5, 2.6 and [T^ for more 
details). 

An isothermal density perturbation of the form ( |l6| ) 
evolves into a regular stream of cooling gas towards the cen- 
tre in the equatorial plane and its close vicinity. The stream 
is embedded into an extended slow outflow region covering 
the rest of the computational volume. The development of 
the instability occurs in several stages. First, the gas den- 
sity starts to grow in the perturbed region near the centre. 
Then, at some point (t ~ 0.5 t s ,o), gas starts to cool and 
condense very efficiently in the inner part of the equatorial 
plane (r < 3 kpc) where the perturbations grow non-linear 
first. This produces a deep and narrow pressure (density) 



minimum sandwiched by sharp pressure (density) maxima 
above and below the equatorial plane, which smoothly ap- 
proach the equilibrium pressure (density) values near the 
angular grid boundaries. The whole configuration works as 
a 'funnel' (hereafter in a 2D sense) which sucks in rarefied 
gas from the periphery of the stratified hydrostatic distribu- 
tion towards the galactic centre. Close to the inner spherical 
boundary the highly accelerated gas stream passes through 
a shock front and quickly condenses out in the dense post- 
shock region. 

In order to model the processes in the hot galactic coro- 
nae more realistically, we superimposed random perturba- 
tions of the velocity field (|18[) onto the equilibrium hydro- 
static gas distribution (P"^)-([14D; see Fig. ^. Since these veloc- 
ity disturbances generate isobaric perturbations (which are 
then being caught by the instability) of considerably smaller 
amplitude, it takes longer to develop a non-linear flow. We 
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^ = 0.000 



£=0.996 



log p 





0.0 0.5 1.0 
log r (kpc) 

Figure 5. Example of an axisymmetric random velocity pertur- 
bation of the form ( |l8| ) with an amplitude e = 10 %. The straight 
line marked with z indicates the symmetry axis. 



find that the flow enters the non-linear regime at t ~ t s 



1.0 t s 



5.0 x 10 yr are presented in 



The results at t 
Figs |0. 

The pressure map (Fig. H) shows three funnels in the 
central region (r < 7 kpc). The number of developed funnels 
and their distribution in angular direction directly reflects 
the number of harmonics present in the initial velocity per- 
turbation and their phases. Non-linear coupling of pertur- 
bation modes gives rise to a different flow development in 
the individual funnels at t = t s ,o- In clockwise order, weak, 
strong, and medium streams can be seen. 

The arrows superimposed on the map give only a rough 
representation of the flow pattern, because the "arrow grid" 
is sparser than the computational grid to avoid illegible fig- 
ures. The velocity field for the central 1 kpc region is shown 
in more detail in Fig. [l(j. Three streams of hot gas flow into 
the low pressure regions located inside the funnels. Simulta- 
neously, some gas leaves the high pressure regions associated 
with the 'funnel walls' and enters the slow outflow regime 
between the funnels. 

The pressure distribution in the leftmost funnel is rather 
patchy. The corresponding features can also be seen in the 
velocity map (Fig. [lo|). The gas tends to concentrate in ra- 
dially compressed clumps, where it condenses out digging a 
low pressure channel into the hot galactic atmosphere. As a 
result a radial hollow jet of high velocity gas forms. Growing 
from linear values of about 100 km s _1 at t = 0.85 t s ,o the 
maximum flow velocity saturates at a slightly supersonic [in 
relation to c(T eq )] value of < 700 km s^ 1 when the temper- 
ature in the funnels drops to ~ 10 6 K. 

The inflow-outflow boundary represents a tangential 
discontinuity. In an ideal fluid such discontinuities are abso- 
lutely u nstable even with respect to infinitely small pertur- 
bations (Landau & Lifshitz 1959). The growth rate of this 



shear instability, lm(u>), depends on the wavenumber of the 
instability, the relative velocity v r of the fluids and on the 
densities of both fluids 



Im(tj) = kv, 



VP1P2 
pi + p-2 ' 




log t (kpc) 

Figure 6. Example of an axisymmetric hybrid inflow-outflow so- 
lution. Displayed are the distribution of gas pressure (colour map) 
and the velocity field (arrows) at 4.9 X 10 7 yr, which forms in re- 
sponse to the random velocity perturbation shown in Fig. 5. The 
colour bar labels give p on a logarithmic scale in arbitrary units. 



Z = 0.996 




T (K) 

n 
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Figure 7. Same as Fig. 6, but showing the gas temperature (on 
a linear scale). 



On time-scales of ~ O.H s ,o perturbations with a wavelength 



A < A cr it = 0.6 



8 P 



lOOkms -1 1 + Sp 



kpc, 



(30) 



(29) 



where Sp = become non-linear. However, our case is 
considerably more complex, since in the non-linear regime 
compressibility effects are important and since the instabil- 
ity of the tangential discontinuity can be coupled with the 
condensational instability driven by the sources. Neverthe- 
less, one can still obtain a rough estimate from equation 
(pp|). For a typical density contrast Sp — 0.1 and a relative 
flow speed v v = 400 km s _1 one finds A cr i t = 0.7 kpc. This 
implies that the instability may be important for the further 
evolution of the streams. 

The temperature map (Fig. shows that the initial 
isothermal gas distribution gets destroyed by the perturba- 
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Figure 8. Same as Fig. 6, but showing the logarithm of the gas 
number density in units of 0.28 cm -3 . 
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Figure 10. Velocity field in the central region at t = 4.9 X 10 7 yr. 
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Figure 9. Same as Fig. 6, but showing the logarithm of the 
normalized deposited gas density. 



tions. The funnels are filled by cooler gas with temperatures 
as low as 10 6 K, while in the outflow regions the gas is hotter 
than in the equilibrium (T « 2x 10 7 K). Note, that although 
periodic boundary conditions are applied in angular direc- 
tion, there exists no high temperature region between the 
leftmost and (the periodically continued) rightmost funnel. 
Instead, there are indications of a temperature drop and in- 
flow in the inter-funnel region. This may be considered as 
an indication of an imminent coalescence of the two evolved 
funnels, because of the same reason as in the case of the coa- 



lescence of condensation fronts described in (Kritsuk 1998). 
However, in the realistic three-dimensional case the proba- 
bility for such a coalescence of two finger-like streams must 
be quite low. It seems that at least the hot, well structured 
outflow region between the central and the rightmost funnel 
will prevent them from merging. 

The (number) density map (Fig. ^) looks qualitatively 
similar to the pressure map. Different stages of gas com- 
pression and deposition can be seen in the three funnels. 



The maximum number density is 1.7 cm . Such large val- 
ues are reached in the condensing fragments falling down 
towards the centre.^] 

The map of the normalized mass density (Fig. ^) de- 
posited during the simulation is similar to the temperature 
map. Gas condensation was much more efficient in the fun- 
nels, while in the outflow regions the rate falls below the 
equilibrium values. The condensation pattern in the form- 
ing streams looks rather patchy. This indicates once again 
that during the development of the streams the gas cools 
and condenses out in clumps. With an initial number den- 
sity of 1.0 cm" 3 and a size of 0.2 kpc the mass of such a 
clump, ~ 10 6 Mq, corresponds to the mass of a globular 
cluster. 

Apparently, the presence of hollow jets, which fill only 
a small fraction of the whole volume in the central parts of a 
galaxy, does not change the X-ray observables dramatically 
as compared to those of the equilibrium gas configuration, if 
individual streams cannot be resolved well. To demonstrate 
this we show in Fig. |ll] angular averages of the same vari- 
ables, which were plotted for spherically symmetric flows in 
Figs hlM The mean temperature is averaged with a weight 
p 2 , while we use a p- weighted mean for the radial velocity. 
We find that, at least at t — 1.0 t s ,o, the density and sur- 
face brightness profiles do not differ substantially from the 
equilibrium ones. The mean temperature drops by ~ 50 % 
towards the centre showing some oscillations due to indi- 
vidual cooling clumps. Weighted mean inflow velocities up 
to ~ 10 2 km s -1 indicate substantial deviations from hy- 
drostatic equilibrium. [The adiabatic sound velocity at the 
equilibrium temperature is c(T cq ) as 555 km s .] At the 
same time, the bulk (~ 50%) of the condensed gas mass is 
deposited in the very center of the galaxy, r < 2 kpc, and 
the mass deposition rate within the cooling radius is com- 
parable to that of the spherically symmetric inflows shown 
in Figs and |. 



" Note, that this is the density of the hot, pressure supported gas 
phase rather than the density in deposited condensates, which is 
much higher. 
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Figure 11. Example of an axisymmetric hybrid inflow-outflow solution at t = 4.9 X 10 7 yr. Displayed are various angular-averaged 
variables. See also Fig. 1. 



These results suggest that in the absence of direct 
gas velocity measurements the interpretation of the surface 
brightness and the temperature profile in terms of the quasi- 
hydrostatic steady-state approximation built into the 'stan- 
dard' cooling flow model can easily result in an erroneous 
spatial distribution of the mass deposition rate M(r). Viola- 
tion of hydrostatic equilibrium in a small fraction of the vol- 
ume, filled with accelerated gas streams, can lead to a con- 
siderable mass flux of gas towards the centre, while the bulk 
of the hot galactic corona remains in a quasi-equilibrium 
state. 

Because of the extremely short cooling time scale in 
the blobs near the centre, we cannot simulate the develop- 
ment of the instability into the non-linear regime for a long 
time. The final fate of the evolving perturbations, thus, re- 
mains unclear at this stage of the study. It is hard to predict 
whether the formation of funnels represents just the onset 



of a forthcoming global cooling catastrophe at the centre 
or whether the funnels are transient features, which will be 
destroyed by hydrodynamic instabilities. Feedback heating 
could play an important role in saturating the instability 
in the non-linear regime, but the definition of feedback is 
uncertain. 



6 DISCUSSION 

We found solutions, which can be classified in four different 
types: 

(i) A subsonic outflow with temperature enhancement in 
the central region ~ 3 with typical gas velocities of 
about 200 km s _1 and a low X-ray luminosity; 

(ii) A quasi-hydrostatic gas configuration with outflow ve- 
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locities ~ 10 km s _ , a nearly isothermal gas distribution, 
and X-ray observables close to those of the recycling model; 

(iii) An unstable spherically symmetric inflow with a X-ray 
luminosity considerably exceeding observed values; 

(iv) A hybrid inflow-outflow solution with a temperature 
depression in the center, mean inflow velocities of ~ 
100 km s~ , and a X-ray brightness distribution similar to 
the one of the recycling model. 

The prototypes are given in Figs ^, ^, ^f^[ and [H], respec- 
tively. 

Since the spherically symmetric non-linear regime of 
catastrophic cooling is unstable, we consider solutions of 
type (iii) as being not realized in nature. The other three 
types represent flow regimes, which can be expected to oc- 
cur in galaxies, [JJ] depending on the degree of suppression of 
thermal conductivity in the hot ISM and on the source of 
the large-scale perturbations. 

According to a simple scenario, the hot ISM evolves 
through several stages with a galaxy free of any gas in the 
beginning: (1) a supersonic galactic wind driven by early 
SN activity, (2) a slow outflow regime, and (3) an i nflow, 
trigger ed b y the central cooling catastrophe; see (Ciotti 
et al. 1991). The evolution is driven by coordinated secu- 
lar changes of the control parameters which determine the 
mass loss rate and the characteristic temperature of the dis- 
tributed heat source in the galaxy, a(t) and To(t). Since the 
study of Ciotti et al. (1991) does not include the mass depo- 
sition due to local thermal instabilities in the ISM, one can 
speculate, that their inherently unstable (either to wind or 
to inflow) outflow regime, probably, corresponds to our, in 
the same sense unstable, equilibrium recycling model. In our 
case the hot ISM in elliptical galaxies could be accumulating 
to fill the potential well and to approach 'thermodynamic' 
equilibrium between sinks and sources, while the SN activity 
and stellar mass loss decrease. The evolution in the outflow 
phase near equilibrium could last for a long time before it 
enters the inflow phase and ends up with a cooling catas- 
trophe. It is expected that this transition can also be trig- 
gered earlier by external or internal processes which disturb 
the ISM of the galaxy. However, a thorough study of the 
stability for evolutionary ISM models with time-dependent 
sources is needed in order to draw rigorous conclusions for 
this non-linear problem. 

Our hybrid inflow-outflow solutions indicate that the 
central cooling catastrophe can develop in the form of ra- 
dially oriented streams of low density peripheral gas pen- 
etrating the subsonic outflow. The vertical orientation of 
(nearly) isolated giant ellipticals is defined by the condition 
£g <§C tc,o- When the thermal and gravitational time scales 
are comparable in the central region, or when to > t c ,a, 
the condensation fronts will be orientated in a random way, 
since gravit y does not focu s the flows. This case is described 
in detail in (Kritsuk 1998). The matter, deposited in narrow 
streams, can show up in optical and UV emission lines as fil- 
amentary nebulae residing within a radius of a few kpc from 
the nucleus. The morphology of such filaments depends both 



II Note, that among these types there is no stable spherically 
symmetric 'cooling flow' realization with a subsonic inflow at all 
radii. 



on the ratio of the characteristic thermal and free fall time 
scales and on the topological properties of the large-scale 
perturbations triggering the instability. A detailed study of 
the dynamics of the cooling catastrophe requires multidi- 
mensional high resolution simulations, which would allow 
one to follow the development of the hydrodynamic insta- 
bilities further into the non-equilibrium epoch. This work is 
in progress. 

Numerical experiments allow us to identify two modes 
for mass deposition in hot galactic coronae. The first, hydro- 
static continuous mode, corresponds to a distributed mass 
deposition, as it occurs in equilibrium models and in slow or 
subsonic outflows. The second, dynamical mode, is related 
to the catastrophic cooling regime. These modes can be re- 
sponsible for the diffuse blue light emission due to contin- 
uously distributed star formation of low efficiency and due 
to massive young star clusters formed in bursts triggered 
by large-scale perturbations. The amount of mass deposited 
in both ways can be comparable. However, the mean pti(r) 
distribution follows the density of the old stellar population 
in the quiescent mode, while it is much more concentrated 
towards the centre in case of the catastrophic mode. 

Currently our study does not include the effects of an 
additional intra-cluster gas component, which typically has a 
considerably higher temperature corresponding to the clus- 
ter gravitational potential. The intra-cluster medium can 
serve as a thermal energy reservoir for the ISM of the galaxy 
and it can essentially modify the gas dynamics through the 
enhanced pressure at the outer boundary. We will address 
these problems in a subsequent publication. On the other 
hand, slow rotation and flattening of the galactic gravita- 
tional potential do introduce finite perturbations to the equi- 
librium model. Hence, one expects a hybrid solution with in- 
flow along the equatorial plane and with polar outflows. Re- 
cent results of D'Ercole & Ciotti (1997) have clearly demon- 
strated that this is indeed the case. 



7 CONCLUSIONS AND PERSPECTIVE 

Using numerical technique we probed the stability of the hy- 
drostatic equilibrium recycling model for hot gaseous coro- 
nae of giant elliptical galaxies with respect to a variety of 
perturbations. 

In the absence of heat conduction the equilibrium ap- 
pears to be unstable on a time scale of the order of the 
thermal time for the gas at the galactic centre. The physi- 
cal reason for this instability is rooted in the properties of 
the source terms. The energy and mass sinks due to local 
thermal instabilities essentially modify the elastic proper- 
ties of the gas. A slow contraction of a gas element amplifies 
the condensation and cooling in it and, thereby, the mass 
inflow into the compressed region. As a result, the density 
rises, providing further growth of losses and accelerating the 
contraction in a runaway regime. A slow expansion of the 
gas element instead decreases the losses and increases the 
specific heating per unit mass, which further decreases the 
losses and accelerates the element expansion. The expansion 
and heating saturate as the gas temperature approaces the 
characteristic temperature of the heat source. Both types of 
unstable behaviour can be identified in our simulations. 

Spherically symmetric perturbations trigger either gas 
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inflow, which ends up with a cooling catastrophe in the core 
of a galaxy, or a subsonic outflow regime, which efficiently 
removes the gas from the central region. In the first case the 
X-ray surface brightness at the centre of the galaxy consid- 
erably exceeds observed values, while in the second case the 
hot gas is practically invisible in X-rays. 

Thermal conduction is able to maintain a stable equi- 
librium only against low-amplitude perturbations. Spheri- 
cally symmetric, global density perturbations with ampli- 
tudes > 10 % remain unstable. In case of inflows the conduc- 
tive heat flux cannot prevent the cooling catastrophe at the 
centre, but can delay it in time considerably. On the other 
hand, when the thermal conductivity is not suppressed, the 
gas velocities of outflows saturate at much lower values. This 
gives rise to quasi-steady-state solutions with an X-ray lu- 
minosity only slightly lower than that of the equilibrium 
model. 

Using axisymmetric perturbations, we were able to 
study the 2D hydrodynamics of the cooling catastrophe for 
the first time. In axial symmetry the system exhibits an 
qualitatively different behaviour. A set of narrow cooling 
gas streams, flowing towards the centre through a global 
subsonic galactic wind, develops in response to random ve- 
locity perturbations of the equilibrium recycling model. In 
the strongly non-linear regime the characteristic averaged 
X-ray observables of such hybrid gas flows (e.g., the surface 
brightness and the temperature profiles) mimic the charac- 
teristics of the initial equilibrium state, while the equilibrium 
itself appears to be considerably violated. 

These results indicate the need for a three-dimensional 
time-dependent treatment of the problem, which would be 
able to reveal non-linear instability saturation mechanisms. 
At the same time, they cast doubts on using steady-state 
'cooling flow' type models, based on the assumption of quasi- 
hydrostatic equilibrium, for recovering mass deposition rates 
in the central galactic regions. 
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